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Abstract 

We propose an approximate second order maximum entropy (M2) model 
for radiative transfer in slab geometry. The model is based on the ansatz 
of the specific intensity in the form of a ^-distribution. This gives us an 
explicit form in its closure. The closure is very close to that of the maxi¬ 
mum entropy, thus an approximation of the M2 model. We prove that the 
new model is globally hyperbolic, sharing most of the advantages of the 
maximum entropy closure. Numerical examples illustrate that it provides 
solutions with satisfactory agreement with the M2 model. 
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1 Introduction 

The radiative transfer equation describes the density of a system of particles inter¬ 
acting with a background medium. It has been widely used in various applications 
such as atmospheric modeling, nuclear engineering and medical imaging. As the 
radiative transfer equation is a problem in very high dimension, deriving a low 
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dimensional model is the first step before further numerieal studies. Most models 
for radiative transfer, and more broadly, for kinetie equations in general, fall in 
one of the following eatogories: partiele models, moment models, and diserete- 
veloeity models. In this paper we will foeus only on moment models. 

Often, the moment models are equipped with the niee property of being natu¬ 
rally rotational invariant, their variables have elear physieal meaning and therefore 
offer elear insight into the physies of the problem under eonsideration. They ean 
be highly effieient in many applieations, for example the Euler and Navier-Stokes 
equations in fluid dynamies. However, when ereating moment model it ean be 
diffleult to ensure that it is both hyperbolie and eonsistent with the faet that the 
unknown in the radiative transfer equation is a density and thus must be nonneg¬ 
ative. This seeond property we eall positivity and is a major problem with the 
method [|7||, a linear method whieh is one of the most well-known moment meth¬ 
ods in the radiative transfer eommunity. It was proven in ifT^ that linear moment 
models that are also hyperbolie and rotational invariant almost inevitably bring 
non-positive solutions, thus motivating the study of nonlinear models. However, 
it is not straightforward to give a hyperbolie nonlinear model0 A hyperbolie mo¬ 
ment model sueeessfully preserving positivity is the maximum entropy model in 
radiative transfer. It was first proposed by Minerbo ifTTI . and Levermore general¬ 
ized it and exposed its mathematieal strueture flUl. Unfortunately, the maximum 
entropy model eurrently has no effieient implementation beeause the elosure rela¬ 
tion is not explieit but must be eomputed through the solution of an optimization 
problem. However, it is still regarded as the most attraetive model due to its highly 
desirable mathematieal properties. 

While some reeent work has been direeted towards developing effieient algo¬ 
rithms for Mn fflU, there has also been inereasing effort devoted to its approxi¬ 
mation. In this work, we investigate a simple ease, in partieular in slab geometry 
and where frequeney dependenee is omitted. This ineludes the single-frequeney 
and gray medium eases. When only the first three moments are speeified, ob¬ 
servations on the speeifie intensity maximizing the Bose-Einstein entropy in both 
oases gives us the expeotation that the speeifie intensity with maximal entropy ean 
be well approximated by a ^-distribution. With suoh a form, an explieit expres¬ 
sion of the fourth moment as a funotion of the first three moments, is obtained, 
resulting in a moment model with explieit elosure. 

Eor this new model, we illustrate that its elosure is very similar to the elosure 

' For instance, the counterpart of model for the Boltzmann equation is Grad’s model im. 
Unlike P„, it is nonlinear, and is not globally hyperbolic. 
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of M 2 model, both for single-frequency and gray case. The three closures are con¬ 
sistent with an underlying nonnegative density and agree exactly on the boundary 
of the realizability region. In the interior of the realizability region, all three clo¬ 
sures agree with each other qualitatively quite well. Moreover, we show that the 
new model shares two important mathematical properties of the M 2 model: global 
hyperbolicity and finite signal speeds no larger than the speed of light. 

The new approximate model is very convenient for numerical simulation due 
to its explicit closure and conservative formulation. We present several numerical 
results comparing with the results given by the original M 2 model. Due to the se¬ 
rious difficulties in the implementation of the M 2 model, the numerical results of 
the M 2 model is obtained by very complex techniques, precisely a revised imple¬ 
mentation based on [|3. Even so, it is still very time consuming. The comparison 
of the numerical results are quite satisfactory, considering the dramatic efficiency 
improvement by the approximate model. 

The rest of this paper is arranged as follows: In Sec. [2]we introduce the basics 
of moment models and the maximum entropy closure. In Sec. |3] we propose 
an approximate M 2 model and analyze its properties. In Sec. |4] we compare 
the approximate model with M 2 using several examples. Finally in Sec. [5] we 
summarize and draw conclusions. 

2 Preliminaries 

Let the specific intensity / (f, u, Ct) to be proportional to the density of radiation 
energy, which is a function depended on the time t E the spatial coordinates 
a; G the frequency u E and the direction variable Ct E on the unit 
sphere. It is governed by the general form of radiative transfer equation as 


1 dl 

+ n.wi = c{i), 

c at 


( 1 ) 


where C{I) describes the interaction of radiation with background medium. 

Denote by {rrii (fi)} a set of basis of a polynomial space Hn ($7), then the 
moments of the specific intensity I are {ml) where we use the notation (■) for 
either 




or 


where dS{Q) is the volume element on the sphere. The former leads only to an 
angular closure, while the latter leads to the so-called gray approximations. The 
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( 2 ) 


exact moments {ml) satsify 

\ f) ITTi T\ 
c ot 

where (g) can either mean g dSl or g dS7) dv. However, the equation 

for the term (flml) involves moments of polynomials not in Hn (f2), therefore 
dH) is not a closed system. A moment model is then defined by approximating the 
higher-order moments in terms of lower order moments to give a closed system of 
equations approximating resulting in a system of equations of the form 

1 

-^ + V-f{E) = r{E), (3) 

c ot 

where E ~ {ml), f{E) ~ {Ctml), and r{E) ~ {mC{I)). How this closure is 
made is called a moment closure and has a fundamental impact on the performance 
of the moment model. 

The maximum entropy principle is an elegant way of deriving moment closure. 
It is based on reconstructing an ansatz of / from the moments by solving the 
following constrained variational maximization problem 

max Hil) 

(4) 

s.t. {Im) = E, 


where H(I) is the Bose-Einstein entropy 




-Jlog(J)+ J + 


2hiy^ 


log I 


2hu^ 


For the angular closure, the solution of dH) has the form 


Lm = 


2hu^ 


hu 

exp I — a ■ m 

Kb 


-1 


while for the gray approximations the solution of dH) has the form 

Lm = "" 


a • m] 


(5) 


( 6 ) 

(7) 


where a is the Stefan-Boltzmann constant. In both cases cx = cx.{E) is the unique 
vector such that {la'm) = E. Then the method is defined by taking 


f{E) = (Q,miJj and r{E) (mC{iS) (8) 
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in ([3]). 

Properties of the model are discussed in [|3 |^, including a proof of its 
global hyperbolicity. It is also positivity preserving, and entropy dissipating. 
However, from dH) one see that the closure is not given explicitly. Instead one 
has to solve (lam) = E for the Lagrange multipliers ol, which involves solving 
a coupled nonlinear algebraic system. Unfortunately, it is expensive and difficult 
to numerically solve this algebraic system. Due to these numerical difficulties, 
there has so far been no efficient general implementation of the Mn model except 
in the Mi case mini. However, in some examples, the Mi model is qualitatively 
wrong [HI . Generally, there are two approaches for resolving the difficulties in 
the implementation of the maximum entropy model. One approach is to develop 
efficient algorithms for solving the optimization problem. There has recently been 
some progress in computing for high order M„ models [[3113- The other approach 
is to give an approximate model of M„ which is explicit and therefore more com¬ 
putationally feasible, while still preserving as many of the advantages of the M„ 
model as possible. 


3 Approximate M 2 Model 

Due to the difficulties in deriving an approximate model, we restrict ourselves to 
the radiative transfer equation in slab geometry and consider only an approxima¬ 
tion of M 2 model. The radiative transfer equation becomes 

/re [-1,1], Io = 


For the case where we only perform an angular closure, /, (Tq, and ag are still 
dependent on the frequency ly, while in the case of the gray approximations, we 
assume that u has already been integrated out of the equation. Therefore from 
now on the angle-bracket notation indicates integrals over /i 



g{fi) d/i. 


Let Ej ~ 
M 2 , then 



/iM(/i) d/i, and denote by Ai the realizable moment vector space of 


Ai — {{Eq, El, E 2 ) |0 < U 2 < Eq, E^ < EqE2 } . 


5 




as shown in IfT^ . We observed that the specific intensity with maximal entropy 
may be well approximated by a ^-distribution if we consider a M 2 model. This 
makes us take the ^-distribution as an ansatz for / : 


f_ -^0 <r_“ 

2 ) \ 2 J ’ 



( 10 ) 


We note that the combinations of .B-distribution was used as an approximation 
for the specific intensity in ifTSll . though for a different purpose. 


Remark 3.1. The choice of B-distribution as ansatz is somewhat arbitrary, but it 
has much flexibility, allowing for skewness and non-symmetry. We will show later 
on that it captures the essential profile of the specific intensity. 


Self-consistency for the first to the third moment require 


EflEo + 1 „ {E^/Eof - E2/E0 

Oi = -, p = 


which gives a non-linear closure as 

^3 = J p^I{p)dp = 


E 2 I Eq — 1 


EflEl + 2 El - 3E0E2 
2El - E 0 E 2 - El 


On the boundaries of M., since there is only one nonnegative ansatz with the 
correct moments |l5l and the B distribution is clearly nonnegative, our closure 
agrees with the M 2 closure. Indeed, 

1. lfE2/Eo = 1, then the specific intensity ansatz of M 2 is / = - {Eq + Ef) S{p— 

1) -I- ^ {Eq — El) 6{p 1), for which B-closure shares with M 2 the same 
closure, and E^ = Ei. 


2. If E 2 /Eq = (Bi/Bq)^, then the specific intensity of M 2 is / = Eo6{p — 
Ei/Eq), for which B-closure also shares the same closure with M 2 , thus 
Ez = EH El 

Furthermore, this closure is correct in the isotropic case (that is, when E = 
{Eq, 0, Bo/3) contains the moments of the constant density I{p) = Bo/2). 

In Figure [H we plot the contours of E^/Eq on Ai for a comparison between 
the M 2 model and the B-closure model. Clearly the models agree qualitatively 
quite well. 
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Figure 1: Comparing contour between M 2 for single frequency (left), M 2 for gray 
approximation (center), and S-closure (right) 


Eq= 1,E,=0 ,£2=1/9 Eq= 1,E^=1/2,£3=2/3 




Figure 2: Comparing specific intensity between M 2 and S-closure model 
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For four sets of moments in typieal regions, Figure [2] eompares the speeifie 
intensity between the M 2 model and our B-elosure model. They all qualitatively 
agree with eaeh other. 

Let us show firstly that the new model based on the B-distribution is globally 
hyperbolie. The Jaeobian matrix of the approximate model is 


/ 


J = 


0 

0 

dEs 


1 

0 

dE, 


0 

1 

dEs 


\ 


( 11 ) 


V dEo dEi dE2 / 


T . ^^3 ^^3 . dEs . . . 

Let oo = 7 :;—, Oi = 7 -—, and 02 = 7 r—, they satisfy 

oEi 


0 


Oq — 


(l\ — 


CL 2 — 


dEi dE2 

El [Eq — E2) (3 EqE 2 — 4 Ei^ + i? 2 ^) 

{Eo^ + EoE2-2Ei^Y 
(3 Eq^ — EqE2 — 2 Ei^^ (^EqE2 — 2 Ei^ + i?2^) 
{Eo^ + EoE2-2Ei^Y 
El [Eq — E2) (3 Eq'^ + EqE2 — 4 


{Eq'^ + EqE2 


2\2 


2Ei‘^) 


The eharaeteristie polynomial of J is 


p(A) — — 02 A — fliA — oq. 


( 12 ) 


We then have the following theorem: 


Theorem 3.1 (Global striet hyperbolieity of S-elosure model). The B-closure 
model is globally strictly hyperbolic in the interior of the realizable region M, 
and its propagation speed is less than the speed of light. 

Proof. Let us study it in two eases: 


1. For El = 0, 


p[\) = A 



E2[3Eq — E 2 ) A 
Eq[Eq + E 2 ) J 


Clearly, the three distinet roots of p[X) are Aq = 0 and X± = 
E2{3Eo — E 2 ) 


Since 0 < 


Eq[Eq + E 2 ) 


< 1, all the roots are within (—1,1). 


E2[3Eq — E 2 ) 
Eq[Eq + E 2 ) 
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2. Consider Ei ^ 0. Without loss of generality, assume Ei > 0. Then 


Oo -El 

- 1 < 0 < — <^< 1 . 
3 Eo 


LaQ = ^- 


EA E 2 

-),Z=l--,.hen 


02 


4 El 


p{^) = ^^Q zf{Q,z)>o, 


with 


f{Q,z) = 


QSQ'^Z^ + 86QZ^ + 27^4 + 288g3 + SGOQ^Z + 112QZ2 

{2QTW 


We notiee that 


_ ^ {Eo — E2 Y{Eq + Ei){Eo + 2Ei + E 2 ) 

’ ~ {Ei + E0E2 - 2 EIY 

(El\ = 4^1 {Eq - Ei){Eo + Ei){EoE 2 - Elf 

^\Eo) El {Ei + EoE2-2E‘lY 

^ {Eo - E2 )\Eo - Ei){Eo - 2Ei + E 2 ) ^ ^ 

=- (Ei + E„E,-2E?r- -> «’ 


thus p{\) has one root in each intervals 
Similar arguments work for i?i < 0. 



/02 EA 
['s^WoJ 


, and 



This ends the proof. □ 

Denote Ai < A 2 < A 3 to be the three eigenvalues of the Jacobian matrix J, 
then it is clear that the eigenvector corresponding to Xj is = [l Aj A^] 

Theorem 3.2. The -characteristic fields are genuinely non-linear, while 

the -characteristic field is neither genuinely non-linear nor linearly degener¬ 

ate^ 

^For the definition of genuiely non-linear and linearly degenerate characteristic fields, see m. 
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Proof. Let 


As 


and 


A(A) = [1 A A^] ■ 


5 (ao, Oi, 02) 
d{Eo, El, E2) 


1 

A 

A 2 


VAj- ■ 




A(A), 


dp 

dX 


> 0, j 

A=A,' 


1,3; 


dp 


< 0 , 

A=A2 


vXj ■ = 0 is equivalent to Xj being the eommon root of A (A) and p{X). As 

the resultant of p{X) and A(A) is 


res(p, A, A) 


^ EiQ^z^ (z2 +16 g + 8 z) (g + z) (z^ + 4 g)" 

izo -Tc- 

E^{2Q + Z)^^ 


E2 (El A E2 

where Q = — -^ , Z = 1 — —. As Q and Z are both positive in the 

Eq \EoJ Eq 

interior of the realizable region M., we have that res(p, A, A) 7 ^ 0 if Zi 7 ^ 0 . 
Therefore when Ei f 0, all eharaeteristie fields satisfy VAj • 7 ^ 0. 

In ease that Zi = 0, A 2 = 0 is a eommon root of p{X) and A(A), so 

VA 2 • = 0 . 


Meanwhile we have 


A(A) = 2A(A2Zo-Z2) 


(3 Eq + E 2 ) {Eq — E 2 ) 
Eq^ {Eq + E2f ’ 


thus VAj ■ R^^^ 7 ^ 0 for j = 1, 3. 

Colleeting the arguments above, one has that 

VA, ■ ^ 0, V(Zo, Ei,E 2 )eM, for j = 1 , 3; 

sign (VAj ■ j = sign(£'i), under proper sealing of R^‘^\ 


□ 


10 










We point out that it is also valid for the M 2 model that VA 2 ■ = 0 when 

El = 0. Aetually, the specific intensity of the M 2 model is 




exp(l + aijj, + a 2 /i^) — 1 ’ 


for the single-frequency case. In this formation, the denominator is always posi¬ 
tive on /i G [—1,1], and we have that Ei = 0 implies ai = 0. Therefore, 


E^ — 0, A 2 — 0, V {Eq, 0, E 2 ) G M.. 



Direct calculations show VA 2 • 1 (Eo,o,E2 ) 


OEq ^ 


Let us summarize briefly some advantages of the new model: 

1. The ansatz for the specific intensity preserves positivity and has explicit 
closure relationship; 

2. The model derived is conservative and globally hyperbolic; 

3. The signal speed of the new model is less than the speed of light; 

4. The closure is very close to that of M 2 ; 

We present some numerical results for some benchmark problems to show the 
quality of the new model as an approximation of the M 2 model. 

4 Numerical Results 

Our approximate M 2 model of dH) is therefore 



(13) 
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We consider the angular closure for a single frequency u = 1. We solve equation 
(fT3]) using the canonical finite volume scheme with the Lax-Friedrich numerical 
flux, and the source term is treated implicitly. 

To impose an inflow boundary condition, we only need to impose the value of 
the flux on the boundary. We derive it using upwind on the kinetic scale. As we 
know the moments on the left and right cells, we can reconstruct I{Eq^ E\^ E 2 , jj) 
and I{Eq, El, E^, /r) on the left and right cells using the ansatz (fTOl) . then integrate 
over to have 

E, = El E[, E[) d/i + J(/r, El E^) d/i, for j = 1, 2, 3, 

to give the flux on the boundary. 

We compare our model to numerical solutions of the true M 2 model, which for 
the single-frequency case uses the ansatz ® . We compute M 2 solutions using the 
kinetic scheme and optimization techniques given in |l2l . The entropy from that 
work is replaced with the Bose-Einstein entropy ([5]), and to avoid the singularity 
in the ansatz ® when the polynomial a. ■ m passes through zero, we limit the 
step-size in the Armijo line search so that the polynomial ct ■ m remains negative 
at every angular quadrature point. Our solutions are computed with 1000 cells, 
and we note that none of the computations below required the use of the isotropic 
regularization technique. 

Below we give the numerical results for three examples. 

Example 4.1 (Two-beam). The absorption coefficient is (Tq = 2, the scattering 
coefficient as = 0, and the speed of light is taken to be c = 1. The spatial domain 

is z G [0,1]. The initial value is set as Eq = 10“^, Ei = 0, E 2 = -E^for all 

o 

2 ;. Inflow boundary condition are imposed on both ends, thus the specific intensity 
on the boundaries are I(t, 0, p) = exp(—10(/i — 1)^) on the left boundary z = 0, 
and I{t, 1,/i) = exp(—10(/i -f 1)^) on the right boundary z = 1. 

In Figure[3]we compares the steady-state solution of Eq between the M 2 model 
and the .B-closure model. 

Example 4.2 (Isotropic inflow into vacuum). In this example, we consider the 
spatial domain z G [—cxo, 1] with an isotropic inflow source is imposed on the right 
boundary into a domain which is unbounded on the left. We take aa = as = 0. 

Initially, for all z, we take Eq = 10“®, Ei = 0, E 2 = -Eq. The isotropic inflow is 

3 

specified at z = 1. The specific intensity outside the right boundary is I(p) = 0.5. 
We carry out the computation from Iq = D tot = 0.5418 and 0.8. 
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0.8 


Figure 3: steady state solution for two-beam problem 

The results are in Figured whieh are the value of Eoatt = 0.5418 and t = 0.8 
for both the S-closure model and the M 2 model. 

Example 4.3 (Plane source). In this test the spatial domain is unbounded, and 

the initial value is taken as Eq{z) = S{z) + 10“^ Ei = 0 and E 2 = -Eq. The 

o 

simulation time interval is from to = 0 to t = 0.7005 and t = 0.9043. 

The numerical results of Eq for M 2 and the S-closure model are in Figure [H 


5 Conclusion 

An approximate M 2 model for the radiative transfer slab geometry in the cases of 
single-frequency and grey medium is proposed. The new model is based on an 
ansatz formulated as a B-distribution. It shares most of the advantages of the M 2 
model while it has an explicit closure. We are now working the extension of this 
idea to a three dimensional configuration and a many moment model. 
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Figure 4: Results for inflow into vaoumm problem 


1=0.7005 



1=0.9043 



Figure 5: Results for plane souree problem 
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